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Abstract 

We address the issue of interaction between zero-point vacancies in solid ^He as described within 
the shadow wave-function model. Applying the reversible-work method and taking into account 
finite-size effects, we obtain a zero-point monovacancy concentration of (2.03 it 0.02)10"^, which 
is slightly higher than the result due to Reatto et al. for the same model. Utilizing the same 
methodology, we then consider the divacancy, taking into account both the in-plane as well as 
out-of-plane configurations with respect to the basal plane. We find no significant anisotropy 
between both conformation. Furthermore, although there is a small binding tendency, the expected 
divacancy concentration is only 4 — 5 times larger than the value expected in the absence of any 
clustering propensity, 2.5 x 10~^. This result suggests that, within the employed model description, 
no vacancy aggregation leading to phase separation is to be expected in the ground state. 

PACS numbers: 67.80.B-,61.72.jd,61.72.Bb 
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The experimental observation of a non-classical rotational inertia (NCRI) by Kim and 
Chanii^ in solid ^He and its interpretation as a possible manifestation of mass superfiow 
in the solid phase, has triggered intensive research efforts both at experimental and theo- 
retical levels^i^i^. While there is still no consensus as to the precise origin of the observed 
phenomenology, it is a generally accepted notion that it is not an intrinsic property of the 
pristine, defect-free crystalline phase, but rather that it is somehow related to crystal disor- 
der, in the form of crystal defects such as vacancies and interstitials, dislocations and grain 
boundaries^'^iSi'^i^iiiii^ii^'i^ and/or the presence of glassy regionsi^ii^. 

In this context, the role of vacancies has received a significant amount of attention ever 
since the first theoretical proposals of a supersolid ground state. Yet, their role remains 
controversial to this day. On the one hand, finite-temperature path-integral Monte Carlo 
(PIMC) calculations seem to suggest that the ground state of solid ^He does not contain 
vacancies^ii^. On the other hand it has been argued that a number of technical issues, 
including the use of periodic boundary conditions and small numbers of particles in the 
simulation box, may in fact prevent such calculations from correctly assessing the true 
nature of the zero-temperature ground state^*^. Moreover, experimental data indicating 
that a vacancy concentration below 0.4% cannot be ruled out,— as well as arguments due 
to Anderson and co-workers^*^, contend the possibility of zero-point vacancies in solid '^He. 
In this light the question whether or not solid ^He contains a finite zero-point vacancy 
concentration remains open. 

In case of the existence of a finite zero-point vacancy concentration, an issue of concern 
involves the interaction between vacancies^*^ and a possible propensity toward clustering. 
Two recent studies^"^ have addressed this issue from two different points of view. Mahan 
and Shin performed a theoretical study of the interaction between two fixed vacancies in solid 
'^He using elasticity theory of the hep crystal. Their results indicate that the interaction is 
attractive and anisotropic: the divacancy interaction energy was found to be more strongly 
attractive along the c-axis compared to directions within the basal plane. These calculations, 
however, do not explicitly include the effects associated with the zero-point motion of the 
vacancies. Rossi and coworkersi^ addressed this issue by performing calculations based on 
the shadow wave function (SWF) description^^i^^ which has bee n able to reproduce many 
of the properties of ^He in the solid phasesM-^. Their approach is based on the well-known 
equivalence between the calculation of the zero-point concentration of point defects in the 
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ground state of a bosonic quantum system and that of an associated classical solid^*^!^ 
at a finite temperature. Employing this approach and the thermodynamic integration (TI) 
technique^i, in which the formation free energy of the vacancy is determined by subtracting 
the free energies of the computational cells with and without vacancy, they determine the 
zero-point monovacancy concentration within the SWF model. The influence of periodic- 
image effects, however, due to the elastic interaction between the periodic images of the 
vacancy, were not taken into account. In addition to the monovacancy, vacancy- vacancy 
interactions were studied by observing systems containing two and three vacan cies and 
collecting statistics in terms of a vacancy-vacancy correlation function. However, the actual 
zero-point divacancy concentration was not determined. 

In this Brief Report we refine the studies carried out by Rossi and coworkers for the 
SWF model of solid ^He in three ways. First, we provide a more accurate result for the 
zero-point monovacancy concentration, taking into account elastic interactions between pe- 
riodic images. To this end we use the reversible-work (RW) method^^, which allows a direct 
computation of the work required to reversibly introduce a vacancy in an initially defect-free 
crystal without the need for the subtraction of large numbers, and a finite-size extrapola- 
tion. This approach was recently applied in the context of bosonic quantum crystals to 
determine the zero-point vacancy concentration of a system described by the Jastrow wave 
function^^. Secondly, applying the same computational scheme, we determine the zero-point 
concentration of the divacancy for the SWF model. Finally, to detect a possible anisotropy 



as in Ref . 
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we distinguish between the in-plane and out-of-plane divacancy configurations 
with respect to the basal plane. 

Formally, the true wave function of a system of bosons in its ground state can be 



written as^^ 



"^{R) = exp 



/qT . (1) 



where R = {ri,r2, . . . ,r7v} stands for the particle coordinates, $ is an effective potential, 
and Qn is a normalization constant. Since ip{R) is positive everywhere, \ip{R)\'^ can, without 
loss of generality, be interpreted as a Boltzmann factor of a classical system described by 
the potential function $(-R) at a temperature ksT = 1. In the SWF variational theory of 
^He,— the ground-state wave function is given as 

'^iR) = ^^jiR) I dS\[e{v,-s.;)iJs{S), (2) 



where the ipj and i/js are the Jastrow factors 



(3) 



with Tii 



Ti — Tjl, and 



N 



ipsiS) = exp - ^ SV{aSij) 



(4) 



i<j 



The latter depends on a rescaled atomic potential V, here taken to be the Aziz poten- 
tial,— and on the distance Sij = \si — Sj\ between auxiliary variables i and j of the 
set S = {si, S2, . . . , stv}. This set of variables is coupled to the particles coordinates R 
through a Gaussian factor, ^(r^ — s^) = exp [— C|rj — Sjp]. The volume integration in 
dS = dsids2 ■ ■ ■ dspi is over the whole space. The variational parameters b, C, 6 and a 
are those that minimize the expectation value of the energy.— 

Within this formulation, the quantum- mechanical probability-density function 
can be associated with the Boltzmann factor of the classical system described by the effective 
potential 



This fictitious system can be thought as composed of interacting trimers, each one com- 
posed of one actual atom and a pair of coupled shadow degrees of freedom. Given the 
aforementioned equivalence between the quantum system and this classical system, the zero- 
point vacancy and divacancy concentrations in the quantum system described by Eq. ([2]) is 
then equal to the thermal equilibrium vacancy and divacancy concentrations in the system 
defined by Eq. ([5]) at a temperature ksT = 1. 

In order to determine these concentrations we now follow the RW method,—!^ which al- 
lows the computation of the formation free energies of the respective defect configurations in 
the fictitious classical system. The RW approach is based on the construction of continuous 
thermodynamic paths that connect a system of interest to a certain reference. In practice 
this is achieved by introducing one or more coupling parameters that measure the progress 
along the given path. By measuring the reversible work along this path one can obtain the 




TV 




(5) 



i<j 
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free-energy difference between the system of interest and the reference. In the particular 
cases of the thermal equilibrium vacancy and divacancy concentrations, we are interested in 
the formation free energies 

AF^ = F{N-1)-^^F{N) 

= ^FiN) + [FiN - 1) - FiN)] (6) 

and 

AFd = F{N - 2) - ^^F{N) 

= lF(iV) + [F(iV-2)-F(iV)], (7) 

respectively. Here F{N) represents the free energy of a defect-free crystal containing N 
atoms, F{N—1) is the free energy of a crystal containing A^— 1 atoms and a monovacancy and 
F{N — 2) is the free energy of a crystal containing N — 2 atoms and two vacancies adjacent to 
each other. According to the second lines of Eq. ([6]) and Eq. ([7]) they can be written in terms 
of the free energy per atom of the defect-free system and the free-energy differences between 
a monovacancy, (divacancy), cell containing — 1, (A^ — 2), atoms and the defect-free cell 
with A^ atoms. Using this partition, we construct separate thermodynamic paths that allow 
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and 



26 the reversible work values 



us to compute both contributions. As detailed in Refs. 
along these paths are measured along finite-time nonequilibrium simulations during which 
the coupling parameters are varied dynamically. In order to eliminate the systematic errors 
associated wit h the nonequilibrium nature of these simulations, the switching processes are 
carried out in both directions. 

To determine the free-energy differences in Eqs. Eq. and Eq. ([7]) we employ thermody- 
namic paths that involves a continuous transformation of the defect-free interacting fictitious 
classical system such that one or more of its atoms are decoupled from the remainder of the 
system. At the same time these atoms are transformed into a collection of noninteracting 
classical harmonic oscillators. Given the specific functional form of the effective classical 



5 



potential in case of the SWF model, we define the thermodynamic path 

N / , \ 5 



U (R,S,S';{A}) = 5^A,aJ — 

i<j ^^^^^ 

N 



N ^ N 

+ 5 A^A, [V{as,,) + V{as[^)] + - \] 

i<j i=l 

M\2 . „. K ^(0)|2 , „. I / ^(0)|2 



X (^K,|r, - r^P + K2|s, - r^P + K2|s: - r^Pj (8) 

in which each trimer i in the system is coupled to a switching parameter Aj, that is allowed 
to vary between and 1. When Aj = 1 for all particles i, the system corresponds to 
the fictitious system described by the effective potential Eq. (jSj). Similarly, when Aj = for 
all i the system corresponds to a collection of 3A^ noninteracting harmonic oscillators. In 
this case, in addition to the degrees of freedom of the atom rj, the coordinates Sj and 
describing the auxiliary degrees of freedom are also transformed into harmonic oscillators, 
centered about the lattice sites r.'^''. The spring constants ki and K2 are chosen such that 
the mean-square displacement of the harmonic oscillators is approximately equal to that of 
the atoms and auxiliary degrees of freedom in the fully interacting system.—"^ 

Using the general path defined by Eq. ([8]) we can determine both contributions in Eq. ([6]) 
and Eq. ([7]). To compute F{N), we set Aj = A for all i, with A varying between and 1. 
The corresponding thermodynamic path describes a transformation from the defect-free 
interacting system into a collection of 3A^ noninteracting classical harmonic oscillators, for 
which the free energy is known analytically. To compute F{N — 1) — F{N), we fix Aj = 1 for 
i ^ k, and switch A^ between and 1. In this manner, the trimer associated with particle k is 
decoupled from the remainder of the system, turning the defect-free interacting system into 
a system containing a monovacancy at lattice site k plus 3 independent harmonic oscillators. 
In a similar fashion, by fixing Aj = 1 for z 7^ A;, / and varying Xk and A; between and 1, one 
inserts two vacancies into the defect-free interacting system, plus 6 independent harmonic 
oscillators. When the lattice sites associated with atoms k and I are nearest neighbors in 
the lattice, this corresponds to the creation of a divacancy. As usual, to avoid singularities 
in the calculation of the reversible work, the vacancies are not allowed to diffuse during the 
reversible creation process, restricting the motion of their nearest-neighbor particles to their 



6 



respective Wigner-Seitz primitive cells.— For convenience the center of mass of the system 
is held fixed during the switching simulations. The volume during the switching simulations 
is held fixed. In order for the results to become independent of this imposed boundary 
condition, a finite-size extrapolation is required. 

In our calculations we apply the Metropolis algorithm to sample configurations from 
We use orthorhombic simulation cells with a hep structure and numbers of atoms varying 
between 180 and 700 at the melting density p = 0.0294A ^ subject to periodic boundary 
conditions. Before starting each switching process the system is equilibrated for at least 
2 X 10^ Monte Carlo (MC) sweeps where 3A^ random attempts are made to move an atom 
or shadow variable. All switching processes are performed using 1.5 x 10"^ MC sweeps per 
process, which is sufficiently slow to guarantee the regime of linear response. The estimates 
of reversible work are obtained as averages over 120 independent forward and backward 
switching processes. 

The results of the calculations are shown in Fig. [1], which shows the monovacancy and 
divacancy formation free energies as a function of the inverse particle number in the computa- 
tional cell. The monovacancy formation free energy is seen to slightly decrease with increas- 
ing system size, which is a consequence of the elastic image interactions due to the periodic 
boundary conditions. The same occurs for the divacancy formation free energies as depicted 
in Panel b). By plotting the results as a function of 1/A^ and extrapolating to iV ^ oo by 
means of linear regressions, we find estimates for the isolated mono and divacancy formation 
free energies. For the monovacancy we obtain AFm = 6.20 ±0.01, in units of kBT=l, which 
gives rise to an equilibrium concentration of Cm = exp(— AFm) = (2.03 ± 0.02) x 10^^. This 
result is slightly higher than the value Cm = (1-4 ± 0.1) x 10"'^ reported by Rossi et al.—. 
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there is 



The results for the divacancy suggest that, in contrast to the findings in Ref. 
no significant anisotropy with respect to the basal plane. The extrapolated formation free 
energies are AF^ = 10.87 ±0.06 and 10.84 ±0.06 for the in-plane and out-of-plane divacancy 
configurations, respectively. Furthermore, comparison with the extrapolated monovacancy 
formation free energy AF^ shows a positive binding free energy Fb = 2Fy — F2v for the 
divacancy, indicating a driving force toward clustering. However, despite the positive value 
of the binding free energy, the corresponding clustering tendency is found to be rather mod- 
est. The divacancy concentration, given by q = (2;/2)c^ exp(Ffe) where z is the number 
of nearest-neighbor sites in the lattice, is only 4-5 times larger than the divacancy concen- 
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tration in the case of a zero binding energy, 2.5 x 10~^. This suggests that, although two 
vacancies do form a bound state, the binding tendency within the present model is not suf- 
ficiently large to provoke vacancy aggregation leading to large-scale phase separation. This 
is in agreement with the observation of Rossi et al.—. 

In summary, we address the issue of interaction between zero-point vacancies in solid 
'^He as described within the shadow wave-function model. Using the reversible-work method 
taking into account finite-size effects, we consider the for both the in-plane and out-of- 
plane configurations with respect to the basal plane. In addition to finding no significant 
anisotropy between both conformation, the expected divacancy concentration is only ~ 4 — 5 
times larger than the value expected in the absence of any clustering propensity. These 
results suggest that, within the employed model description, no vacancy aggregation leading 
to phase separation is to be expected in the ground state. 

The authors acknowledge financial support from the Brazilian agencies FAPESP, CNPq 
and CAPES. Part of the computations were performed at the CENAPAD high-performance 
computing facility at Universidade Estadual de Campinas. 
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FIG. 1: (Color online) Vacancy formation free energies as a function of inverse particle number in 
the computational cell. Free energies are measured in units of ksT = 1. a) Monovacancy. Symbols 
with error bars represent RW data. Line denotes linear regression, b) Divacancy. Triangles 
represent RW data obtained for divacancy perpendicular to basal plane. Square denote RW data 
for divacancy in the basal plane. Dashed and full lines denotes the respective linear regressions. 
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